/*=================================================================================
 * The functions of Reaction
 *=================================================================================*/

/*---------------------------------------------------------------------------------
 * Standard C++ library headers
 *---------------------------------------------------------------------------------*/

#include <fstream>
#include <iostream>
#include <string>
#include <sstream>
#include <iomanip>
#include <omp.h>

/*---------------------------------------------------------------------------------
 * Local headers
 *---------------------------------------------------------------------------------*/

#include "Reaction.h"
#include "FileReader.h"
#include "Function.h"
#include "Global.h"
// #include "cantera/thermo.h"
// #include "cantera/kinetics.h"
// #include "cantera/base/Solution.h"
// #include "cantera/zeroD/IdealGasReactor.h"
// #include "cantera/zeroD/ReactorNet.h"
// #include "cantera/base/Array.h"
// #include "cantera/zerodim.h"
// #include "cantera/thermo/IdealGasPhase.h"
// #include "cantera/numerics/Integrator.h"

using namespace std;

/*---------------------------------------------------------------------------------
    Name:		Reaction::ReactionConstruction

    Input(0):	None

    Function:	Constuct arrays to store information of the reaction

    Return:		None
 *---------------------------------------------------------------------------------*/

void Reaction::ReactionConstruction()
{
    Stoi_F.Initial(NS, NR);
    Stoi_B.Initial(NS, NR);

    Af.Initial(NR);
    Bf.Initial(NR);
    Eaf.Initial(NR);

    React_TB.Initial(NS, NR);

    Mr.Initial(NS);
    Mi.Initial(NS);
    Yi.Initial(NS);

    Coeff0.Initial(9, NS);
    Coeff1.Initial(9, NS);
    Mw.Initial(NS);
    Ri.Initial(NS);

    Hi.Initial(NS);
    Si.Initial(NS);
    Gi.Initial(NS);

    KF.Initial(NR);
    KB.Initial(NR);
    Kp.Initial(NR);
    Kc.Initial(NR);

    RR_F.Initial(NR);
    RR_B.Initial(NR);
    R_TB.Initial(NR);
    RR.Initial(NR);
    Mc.Initial(ni + 2 * bc, nj + 2 * bc, NS);
    Di.Initial(ni + 2 * bc, nj + 2 * bc, NS);
    Wi.Initial(NS);
    CMS.Initial(ni + 2 * bc, nj + 2 * bc, NS);

    WJH1.Initial(NS, NR);
    WJH2.Initial(NS, NR);
    MD.Initial(ni + 2 * bc, nj + 2 * bc, NS + 3, NS + 3);

    P.Initial(NS);
    Q.Initial(NS);

    Nchem.Initial(ni + 2 * bc, nj + 2 * bc);
}

/*---------------------------------------------------------------------------------
    Name:		ifstream& seek_to_line

    Input(2):	The file; The selected line

    Function:	Locate to the selected line

    Return:		None
 *---------------------------------------------------------------------------------*/

ifstream &seek_to_line(ifstream &in, int line)
{
    int i;
    char buf[1024];

    /* Locate to the beginning line of the file */
    in.seekg(0, ios::beg);

    /* Locate to the selected line */
    for (i = 0; i < line; i++)
    {
        in.getline(buf, sizeof(buf));
    }

    return in;
}

/*---------------------------------------------------------------------------------
    Name:		Reaction::ReactionRead

    Input(4):	The mesh nodes; the reaction model; the thermo file

    Function:	Read the reaction model and thermo file;
                Construct arrays to store the information

    Return:		None
 *---------------------------------------------------------------------------------*/

void Reaction::ReactionRead(char *reaction_model, char *thermofile, Array<double> xnode, Array<double> ynode, int bc)
{
    this->xnode = xnode;
    this->ynode = ynode;
    this->bc = bc;

    ni = xnode.GetSize() - 2 * bc;
    nj = ynode.GetSize() - 2 * bc;

    ifstream fin;
    fin.open(reaction_model);
    int a[100] = {0}, i = 0, line = 2;
    string s, item;
    seek_to_line(fin, line);
    getline(fin, s, '\n');

    /* Transfer string to number */
    stringstream text_stream(s);

    /* Read the number of species, reactions and three-body reactions */
    while (std::getline(text_stream, item, '\t'))
    {
        a[i] = stoi(item);
        i++;
    }
    NS = a[0];
    NR = a[1];
    TB = a[2];
    // cout << NS << '\t' << NR << '\t' << TB << endl;
    // cout << endl;

    /* Construction */
    ReactionConstruction();

    /* Read the stoichiometric coefficient */
    line = line + 6;
    seek_to_line(fin, line);

    for (int j = 0; j < NR; j++)
    {
        i = 0;
        getline(fin, s, '\n');
        stringstream reaction_forward(s);
        while (std::getline(reaction_forward, item, '\t'))
        {
            Stoi_F(i, j) = stoi(item);
            // cout << Stoi_F(i, j) << '\t';
            i++;
        }
        // cout << endl;
        i = 0;
        getline(fin, s, '\n');
        stringstream reaction_backward(s);
        while (std::getline(reaction_backward, item, '\t'))
        {
            Stoi_B(i, j) = stoi(item);
            // cout << Stoi_B(i, j) << '\t';
            i++;
        }
        // cout << endl;
    }
    // cout << endl;

    /* Read the pre-exponent factor, temperature index and activation energy */
    line = line + 2 * NR + 2;
    seek_to_line(fin, line);

    for (int j = 0; j < NR; j++)
    {
        i = 0;
        double b[4] = {0};
        while (fin >> s)
        {
            stringstream geek(s);
            geek >> b[i];
            i++;
            if (i == 4)
                break;
        }
        Af(j) = b[1];
        Bf(j) = b[2];
        Eaf(j) = b[3];
        // cout  << Af(j) << '\t' << Bf(j) << '\t' << Eaf(j) << endl;
        // cout << endl;
    }
    // cout << endl;

    /* Read the coefficient of three-body reactions */
    line = line + NR + 3;
    seek_to_line(fin, line);

    for (int j = 0; j < TB; j++)
    {
        i = 0;
        while (fin >> s)
        {
            stringstream geek(s);
            geek >> a[i];
            i++;
            if (i == NS + 1)
                break;
        }

        for (int k = 0; k < NS; k++)
        {
            React_TB(k, a[0] - 1) = a[k + 1];
            // cout  << React_TB(k, a[0]-1) << '\t';
        }
        // cout << endl;
    }
    // cout << endl;

    /* Read the component mole ratio */
    line = line + TB + 2;
    seek_to_line(fin, line);
    i = 0;
    double c[100] = {0};
    while (fin >> s)
    {
        stringstream geek(s);
        geek >> c[i];
        i++;
        if (i == NS)
            break;
    }
    for (int i = 0; i < NS; i++)
    {
        Mr(i) = c[i];
        cout << Mr(i) << '\t';
    }
    cout << endl;
    fin.close();

    ifstream in;
    in.open(thermofile);

    /* Read the NASA Thermochemical polynomial coefficient and component molecular weight */
    line = 15;
    double aa[100] = {0};

    for (int j = 0; j < NS; j++)
    {
        seek_to_line(in, line);
        i = 0;
        while (in >> s)
        {
            stringstream geek(s);
            geek >> aa[i];
            // cout << aa[i] << endl;
            i++;
            if (i == 19)
                break;
        }
        for (int k = 0; k < 9; k++)
        {
            Coeff0(k, j) = aa[k];
            // cout << setiosflags(ios::scientific) << setprecision(9) <<Coeff0(k, j) << '\t';
        }
        // cout << endl;
        for (int k = 0; k < 9; k++)
        {
            Coeff1(k, j) = aa[k + 9];
            // cout << setiosflags(ios::scientific) << setprecision(9)<< Coeff1(k, j) << '\t';
        }
        // cout << endl;
        Mw(j) = aa[18];
        // cout << setiosflags(ios::scientific) << setprecision(7) << Mw(j) << endl;
        line = line + 6;
    }
    fin.close();

    ReactionInitial();

    return;
}

/*---------------------------------------------------------------------------------
    Name:		Reaction::ReactionInitial

    Input(0):	None

    Function:	Calculate the component gas constant, the component mole fraction
                and the mass fraction

    Return:		None
 *---------------------------------------------------------------------------------*/

void Reaction::ReactionInitial()
{
    Function Fun;
    double Total_Mw = Fun.sum(2, Mr, Mw);

    for (int i = 0; i < NS; i++)
    {
        Ri(i) = R / Mw(i);
        Mi(i) = Mr(i) / (Mr.Sum());
        Yi(i) = Mr(i) * Mw(i) / Total_Mw;
    }
}

/*---------------------------------------------------------------------------------
    Name:		Reaction::Trapezoid

    Input(4):	The component mole concentration; The density fraction;
                The temperature; The time step

    Function:	Simulate the chemical reaction process by Trapezoid formula
                Calculate the component mass production rate and density fraction

    Return:		The density fraction
 *---------------------------------------------------------------------------------*/

Array<double> Reaction::Trapezoid(Array<double> Mc_temp, Array<double> Di_temp, Array<double> T, double dt)
{
    for (int j = bc; j < nj + bc; j++)
        for (int i = bc; i < ni + bc; i++)
        {
            T_local = T(i, j);
            for (int k = 0; k < NS; k++)
            {
                Di(i, j, k) = Di_temp(i, j, k) * 1e-3;
                Mc(i, j, k) = Mc_temp(i, j, k) * 1e-6;
                Hi(k) = GetHi(T_local, R, k, Coeff0, Coeff1);
                Si(k) = GetSi(T_local, R, k, Coeff0, Coeff1);
                Gi(k) = Hi(k) - Si(k) * T_local;
            }

            /* Calculate the net reaction rate */
            RR_F.Fill(1.0);
            RR_B.Fill(1.0);
            R_TB.Fill(0);
            for (int r = 0; r < NR; r++)
            {
                double Xr = 0.0, delta_g = 0.0;
                KF(r) = Af(r) * pow(T_local, Bf(r)) * exp(-Eaf(r) / (Ru * T_local));
                for (int k = 0; k < NS; k++)
                {
                    Xr += (Stoi_B(k, r) - Stoi_F(k, r));
                    delta_g += (Stoi_B(k, r) - Stoi_F(k, r)) * Gi(k);
                    RR_F(r) *= pow(Mc(i, j, k), Stoi_F(k, r));
                    RR_B(r) *= pow(Mc(i, j, k), Stoi_B(k, r));
                    R_TB(r) += Mc(i, j, k) * React_TB(k, r);
                }
                Kp(r) = exp(-delta_g / (R * T_local)) * pow((P0 * 1e-6), Xr);
                Kc(r) = Kp(r) * pow(R * T_local, -Xr);
                KB(r) = KF(r) / Kc(r);
                if (R_TB(r) == 0.0)
                    R_TB(r) = 1.0;
                RR(r) = KF(r) * RR_F(r) - KB(r) * RR_B(r);
            }

            /* Calculate the mass prodcution rate */
            for (int k = 0; k < NS; k++)
            {
                double temp = 0.0;
                for (int r = 0; r < NR; r++)
                    temp += (Stoi_B(k, r) - Stoi_F(k, r)) * R_TB(r) * RR(r);
                Wi(k) = Mw(k) * temp;
                CMS(i, j, k) = Wi(k) * 1e3;
            }

            P.Fill(0), Q.Fill(0);
            for (int k = 0; k < NS; k++)
            {
                for (int r = 0; r < NR; r++)
                {
                    Q(k) += (Stoi_B(k, r) * KF(r) * RR_F(r) + Stoi_F(k, r) * KB(r) * RR_B(r)) * R_TB(r);
                    P(k) += (Stoi_B(k, r) * KB(r) * RR_B(r) + Stoi_F(k, r) * KF(r) * RR_F(r)) * R_TB(r);
                }
                Q(k) = Mw(k) * Q(k);
                if (Mc(i, j, k) == 0)
                {
                    P(k) = 0.0;
                }
                else
                {
                    P(k) = P(k) / Mc(i, j, k);
                }
                Di(i, j, k) = ((1.0 - dt / 2.0 * P(k)) * Di(i, j, k) + dt * Q(k)) / (1.0 + dt / 2.0 * P(k)) * 1e3;
            }
        }

    return Di;
}

/* Adaptive */
Array<double> Reaction::Trapezoid(Array<double> Mc_temp, Array<double> Di_temp, Array<double> Yi_temp, Array<double> T, double dtm, int step)
{
    double dt = 0.0, dt_temp = 0.0;

    for (int j = bc; j < nj + bc; j++)
        for (int i = bc; i < ni + bc; i++)
        {
            if (step == 0 || step < Nchem(i, j))
            {
                dt = dtm;
                if (step > 0)
                    dt = dtm / Nchem(i, j);

                T_local = T(i, j);
                for (int k = 0; k < NS; k++)
                {
                    Di(i, j, k) = Di_temp(i, j, k) * 1e-3;
                    Mc(i, j, k) = Mc_temp(i, j, k) * 1e-6;
                    Hi(k) = GetHi(T_local, R, k, Coeff0, Coeff1);
                    Si(k) = GetSi(T_local, R, k, Coeff0, Coeff1);
                    Gi(k) = Hi(k) - Si(k) * T_local;
                }

                /* Calculate the net reaction rate */
                RR_F.Fill(1.0);
                RR_B.Fill(1.0);
                R_TB.Fill(0);
                for (int r = 0; r < NR; r++)
                {
                    double Xr = 0.0, delta_g = 0.0;
                    KF(r) = Af(r) * pow(T_local, Bf(r)) * exp(-Eaf(r) / (Ru * T_local));
                    for (int k = 0; k < NS; k++)
                    {
                        Xr += (Stoi_B(k, r) - Stoi_F(k, r));
                        delta_g += (Stoi_B(k, r) - Stoi_F(k, r)) * Gi(k);
                        RR_F(r) *= pow(Mc(i, j, k), Stoi_F(k, r));
                        RR_B(r) *= pow(Mc(i, j, k), Stoi_B(k, r));
                        R_TB(r) += Mc(i, j, k) * React_TB(k, r);
                    }
                    Kp(r) = exp(-delta_g / (R * T_local)) * pow((P0 * 1e-6), Xr);
                    Kc(r) = Kp(r) * pow(R * T_local, -Xr);
                    KB(r) = KF(r) / Kc(r);
                    if (R_TB(r) == 0.0)
                        R_TB(r) = 1.0;
                    RR(r) = KF(r) * RR_F(r) - KB(r) * RR_B(r);
                }

                /* Calculate the mass prodcution rate */
                for (int k = 0; k < NS; k++)
                {
                    double temp = 0.0;
                    for (int r = 0; r < NR; r++)
                        temp += (Stoi_B(k, r) - Stoi_F(k, r)) * R_TB(r) * RR(r);
                    Wi(k) = Mw(k) * temp;
                    CMS(i, j, k) = Wi(k) * 1e3;
                }

                /* Calculate dt for each mesh node */
                if (step == 0)
                {
                    dt = dtm;
                    for (int k = 0; k < NS; k++)
                    {
                        if (Yi_temp(i, j, k) >= 1e-6)
                        {
                            if (Wi(k) != 0)
                            {
                                dt_temp = abs(-Di(i, j, k) / Wi(k));
                                if (dt_temp < dt)
                                    dt = dt_temp;
                            }
                        }
                    }
                    Nchem(i, j) = int(ceil(dtm / dt));

                    if (Nchem(i, j) <= 0)
                    {
                        cout << "wrong chem " << i << '\t' << j << '\t' << Nchem(i, j) << '\n';
                        abort();
                    }
                }

                P.Fill(0), Q.Fill(0);
                for (int k = 0; k < NS; k++)
                {
                    for (int r = 0; r < NR; r++)
                    {
                        Q(k) += (Stoi_B(k, r) * KF(r) * RR_F(r) + Stoi_F(k, r) * KB(r) * RR_B(r)) * R_TB(r);
                        P(k) += (Stoi_B(k, r) * KB(r) * RR_B(r) + Stoi_F(k, r) * KF(r) * RR_F(r)) * R_TB(r);
                    }
                    Q(k) = Mw(k) * Q(k);
                    if (Mc(i, j, k) == 0)
                    {
                        P(k) = 0.0;
                    }
                    else
                    {
                        P(k) = P(k) / Mc(i, j, k);
                    }
                    Di(i, j, k) = ((1.0 - dt / 2.0 * P(k)) * Di(i, j, k) + dt * Q(k)) / (1.0 + dt / 2.0 * P(k)) * 1e3;
                }
            }
            if (step == 0)
                NchemNow = Nchem.MaxValue();
        }

    return Di;
}

/*---------------------------------------------------------------------------------
    Name:		Reaction::Diagonalized

    Input(4):	The component mole concentration; The density fraction;
                The temperature; The partial T

    Function:	Simulate the chemical reaction process by IMEX method
                Calculate the component mass production rate and diagonal of the Jacobi matrix

    Return:		None
 *---------------------------------------------------------------------------------*/

void Reaction::Diagonalized(Array<double> Mc_temp, Array<double> Di_temp, Array<double> T, Array<double> Partial_T)
{
    for (int j = bc; j < nj + bc; j++)
        for (int i = bc; i < ni + bc; i++)
        {
            T_local = T(i, j);
            for (int k = 0; k < NS; k++)
            {
                Di(i, j, k) = Di_temp(i, j, k) * 1e-3;
                Mc(i, j, k) = Mc_temp(i, j, k) * 1e-6;
                Hi(k) = GetHi(T_local, R, k, Coeff0, Coeff1);
                Si(k) = GetSi(T_local, R, k, Coeff0, Coeff1);
                Gi(k) = Hi(k) - Si(k) * T_local;
            }

            /* Calculate the net reaction rate */
            RR_F.Fill(1.0);
            RR_B.Fill(1.0);
            R_TB.Fill(0);
            for (int r = 0; r < NR; r++)
            {
                double Xr = 0.0, delta_g = 0.0;
                KF(r) = Af(r) * pow(T_local, Bf(r)) * exp(-Eaf(r) / (Ru * T_local));
                for (int k = 0; k < NS; k++)
                {
                    Xr += (Stoi_B(k, r) - Stoi_F(k, r));
                    delta_g += (Stoi_B(k, r) - Stoi_F(k, r)) * Gi(k);
                    RR_F(r) *= pow(Mc(i, j, k), Stoi_F(k, r));
                    RR_B(r) *= pow(Mc(i, j, k), Stoi_B(k, r));
                    R_TB(r) += Mc(i, j, k) * React_TB(k, r);
                }
                Kp(r) = exp(-delta_g / (R * T_local)) * pow((P0 * 1e-6), Xr);
                Kc(r) = Kp(r) * pow(R * T_local, -Xr);
                KB(r) = KF(r) / Kc(r);
                if (R_TB(r) == 0.0)
                    R_TB(r) = 1.0;
                RR(r) = KF(r) * RR_F(r) - KB(r) * RR_B(r);
            }

            /* Calculate the mass prodcution rate */
            for (int k = 0; k < NS; k++)
            {
                double temp = 0.0;
                for (int r = 0; r < NR; r++)
                    temp += (Stoi_B(k, r) - Stoi_F(k, r)) * R_TB(r) * RR(r);
                Wi(k) = Mw(k) * temp;
                CMS(i, j, k) = Wi(k) * 1e3;
            }
            /* Calculate the diagonal of the Jacobi matrix */
            for (int r = 0; r < NR; r++)
                for (int k = 0; k < NS; k++)
                {
                    if (Di(i, j, k) == 0)
                    {
                        WJH1(k, r) = 0.0;
                    }
                    else
                    {
                        WJH1(k, r) = KF(r) * Stoi_F(k, r) * RR_F(r) / Di(i, j, k) - KB(r) * Stoi_B(k, r) * RR_B(r) / Di(i, j, k);
                    }
                    WJH2(k, r) = RR(r) * (Bf(r) / T_local + Eaf(r) / (Ru * pow(T_local, 2)));
                }

            for (int k = 0; k < NS; k++)
            {
                double temp = 0.0;
                for (int r = 0; r < NR; r++)
                    temp += (Stoi_B(k, r) - Stoi_F(k, r)) * ((WJH1(k, r) + WJH2(k, r) * Partial_T(i, j, k)) * R_TB(r) + RR(r) * React_TB(k, r) / Mw(k));
                MD(i, j, k, k) = Mw(k) * temp;
            }
        }
}

/*---------------------------------------------------------------------------------
    Name:		Reaction::FullMatrix

    Input(4):	The component mole concentration; The density fraction;
                The temperature; The partial T

    Function:	Simulate the chemical reaction process by IMEX method
                Calculate the component mass production rate and the exact Jacobian matrix

    Return:		None
 *---------------------------------------------------------------------------------*/

void Reaction::FullMatrix(Array<double> Mc_temp, Array<double> Di_temp, Array<double> T, Array<double> Partial_T)
{
    for (int j = bc; j < nj + bc; j++)
        for (int i = bc; i < ni + bc; i++)
        {
            T_local = T(i, j);
            for (int k = 0; k < NS; k++)
            {
                Di(i, j, k) = Di_temp(i, j, k) * 1e-3;
                Mc(i, j, k) = Mc_temp(i, j, k) * 1e-6;
                Hi(k) = GetHi(T_local, R, k, Coeff0, Coeff1);
                Si(k) = GetSi(T_local, R, k, Coeff0, Coeff1);
                Gi(k) = Hi(k) - Si(k) * T_local;
            }

            /* Calculate the net reaction rate */
            RR_F.Fill(1.0);
            RR_B.Fill(1.0);
            R_TB.Fill(0);
            for (int r = 0; r < NR; r++)
            {
                double Xr = 0.0, delta_g = 0.0;
                KF(r) = Af(r) * pow(T_local, Bf(r)) * exp(-Eaf(r) / (Ru * T_local));
                for (int k = 0; k < NS; k++)
                {
                    Xr += (Stoi_B(k, r) - Stoi_F(k, r));
                    delta_g += (Stoi_B(k, r) - Stoi_F(k, r)) * Gi(k);
                    RR_F(r) *= pow(Mc(i, j, k), Stoi_F(k, r));
                    RR_B(r) *= pow(Mc(i, j, k), Stoi_B(k, r));
                    R_TB(r) += Mc(i, j, k) * React_TB(k, r);
                }
                Kp(r) = exp(-delta_g / (R * T_local)) * pow((P0 * 1e-6), Xr);
                Kc(r) = Kp(r) * pow(R * T_local, -Xr);
                KB(r) = KF(r) / Kc(r);
                if (R_TB(r) == 0.0)
                    R_TB(r) = 1.0;
                RR(r) = KF(r) * RR_F(r) - KB(r) * RR_B(r);
            }

            /* Calculate the mass prodcution rate */
            for (int k = 0; k < NS; k++)
            {
                double temp = 0.0;
                for (int r = 0; r < NR; r++)
                    temp += (Stoi_B(k, r) - Stoi_F(k, r)) * R_TB(r) * RR(r);
                Wi(k) = Mw(k) * temp;
                CMS(i, j, k) = Wi(k) * 1e3;
            }
            /* Calculate the diagonal of the Jacobi matrix */
            for (int r = 0; r < NR; r++)
                for (int k = 0; k < NS; k++)
                {
                    if (Di(i, j, k) == 0)
                    {
                        WJH1(k, r) = 0.0;
                    }
                    else
                    {
                        WJH1(k, r) = KF(r) * Stoi_F(k, r) * RR_F(r) / Di(i, j, k) - KB(r) * Stoi_B(k, r) * RR_B(r) / Di(i, j, k);
                    }
                    WJH2(k, r) = RR(r) * (Bf(r) / T_local + Eaf(r) / (Ru * pow(T_local, 2)));
                }

            for (int k = 0; k < NS; k++)
            {
                double temp = 0.0;
                for (int r = 0; r < NR; r++)
                    temp += (Stoi_B(k, r) - Stoi_F(k, r)) * ((WJH1(k, r) + WJH2(k, r) * Partial_T(i, j, k)) * R_TB(r) + RR(r) * React_TB(k, r) / Mw(k));
                MD(i, j, k, k) = Mw(k) * temp;
            }
        }
}

/*---------------------------------------------------------------------------------
    Name:		Reaction::ReactionS

    Input(4):	The solver type; The density fraction;
                The temperature; The pressure

    Function:	Simulate the chemical reaction process
                Calculate the component mass production rate
                difference is that the solver is Cantera or DNN

    Return:		The updated component density fraction
 *---------------------------------------------------------------------------------*/

Array<double> Reaction::ReactionS(int type, Array<double> Di_temp, Array<double> T, Array<double> P, double dt)
{
    // double frac_t;
    // double tot_t = dt;
    // double sub_t = 0.0;
    // double T_local, P_local, D_local;

    // double *Di_local = new double[NS];
    // if (type == 0)
    // {
    //     auto sol = Cantera::newSolution("JMH2.yaml", "gas", "none");
    //     auto gas = sol->thermo();
    //     Cantera::IdealGasReactor r;
    //     for (int j = bc; j < nj + bc; j++)
    //     {
    //         for (int i = bc; i < ni + bc; i++)
    //         {
    //             D_local = 0.0;
    //             for (int k = 0; k < NS; k++)
    //             {
    //                 Di_local[k] = Di_temp(i, j, k);
    //                 D_local += Di_local[k];
    //             }
    //             sub_t = 0.0;
    //             T_local = T(i, j);
    //             P_local = P(i, j);
    //             // different frac_t for different zones
    //             if (T_local < 500)
    //                 frac_t = 5e-9;
    //             else
    //                 frac_t = 1e-9;

    //             for (int k = 0; k < NS; k++)
    //             {
    //                 Di_local[k] = Di_local[k] / D_local;
    //             }

    //             gas->setState_TPY(T_local, P_local, Di_local);
    //             r.insert(sol);
    //             Cantera::ReactorNet net;
    //             net.addReactor(r);
    //             net.setMaxTimeStep(frac_t);
    //             net.advance(dt);

    //             for (int s = 0; s < NS; s++)
    //             {
    //                 Di(i, j, s) = r.massFraction(s) * r.density();
    //             }
    //         }
    //     }
    // }
    // else if (type == 1) // DNN
    // {
    // }

    // delete[] Di_local;

    return Di;

}

/*---------------------------------------------------------------------------------
    Name:		Reaction::GetCpi

    Input(3):	Temperature; Gas constant; The numbering of the component

    Function:	Get the specific heat at constant pressure for the given condition

    Return:		The component Cp for the given condition
 *---------------------------------------------------------------------------------*/

double Reaction::GetCpi(double T, double R, int SP, Array<double> Coeff0, Array<double> Coeff1)
{
    double Cpi = 0.0;

    if (T < 1000)
        Cpi = R * (Coeff0(0, SP) * pow(T, -2) + Coeff0(1, SP) * pow(T, -1) + Coeff0(2, SP) + Coeff0(3, SP) * T + Coeff0(4, SP) * pow(T, 2) + Coeff0(5, SP) * pow(T, 3) + Coeff0(6, SP) * pow(T, 4));
    else
        Cpi = R * (Coeff1(0, SP) * pow(T, -2) + Coeff1(1, SP) * pow(T, -1) + Coeff1(2, SP) + Coeff1(3, SP) * T + Coeff1(4, SP) * pow(T, 2) + Coeff1(5, SP) * pow(T, 3) + Coeff1(6, SP) * pow(T, 4));

    return Cpi;
}

/*---------------------------------------------------------------------------------
    Name:		Reaction::GetHi

    Input(3):	Temperature; Gas constant; The numbering of the component

    Function:	Get the enthalpy for the given condition

    Return:	The component enthalpy for the given condition
 *---------------------------------------------------------------------------------*/

double Reaction::GetHi(double T, double R, int SP, Array<double> Coeff0, Array<double> Coeff1)
{
    double Hi = 0.0;

    if (T < 1000)
        Hi = R * (-Coeff0(0, SP) * pow(T, -1) + Coeff0(1, SP) * log(T) + Coeff0(2, SP) * T + Coeff0(3, SP) * pow(T, 2) / 2 + Coeff0(4, SP) * pow(T, 3) / 3 + Coeff0(5, SP) * pow(T, 4) / 4 + Coeff0(6, SP) * pow(T, 5) / 5 + Coeff0(7, SP));
    else
        Hi = R * (-Coeff1(0, SP) * pow(T, -1) + Coeff1(1, SP) * log(T) + Coeff1(2, SP) * T + Coeff1(3, SP) * pow(T, 2) / 2 + Coeff1(4, SP) * pow(T, 3) / 3 + Coeff1(5, SP) * pow(T, 4) / 4 + Coeff1(6, SP) * pow(T, 5) / 5 + Coeff1(7, SP));

    return Hi;
}

/*---------------------------------------------------------------------------------
    Name:		Reaction::GetSi

    Input(3):	Temperature; Gas constant; The numbering of the component

    Function:	Get the entropy for the given condition

    Return:		The component entropy for the given condition
 *---------------------------------------------------------------------------------*/

double Reaction::GetSi(double T, double R, int SP, Array<double> Coeff0, Array<double> Coeff1)
{
    double Si = 0.0;

    if (T < 1000)
        Si = R * (-Coeff0(0, SP) * pow(T, -2) / 2 - Coeff0(1, SP) * pow(T, -1) + Coeff0(2, SP) * log(T) + Coeff0(3, SP) * T + Coeff0(4, SP) * pow(T, 2) / 2 + Coeff0(5, SP) * pow(T, 3) / 3 + Coeff0(6, SP) * pow(T, 4) / 4 + Coeff0(8, SP));
    else
        Si = R * (-Coeff1(0, SP) * pow(T, -2) / 2 - Coeff1(1, SP) * pow(T, -1) + Coeff1(2, SP) * log(T) + Coeff1(3, SP) * T + Coeff1(4, SP) * pow(T, 2) / 2 + Coeff1(5, SP) * pow(T, 3) / 3 + Coeff1(6, SP) * pow(T, 4) / 4 + Coeff1(8, SP));

    return Si;
}

/*---------------------------------------------------------------------------------
    Name:		Reaction::GetWi

    Input(0):	None

    Function:	Get the component mass production rate

    Return:		The component mass production rate
 *---------------------------------------------------------------------------------*/

Array<double> Reaction::GetWi()
{
    return CMS;
}

/*---------------------------------------------------------------------------------
    Name:		Reaction::GetMD

    Input(0):	None

    Function:	Get the diagonal of the Jacobi matrix

    Return:		The diagonal of the Jacobi matrix
 *---------------------------------------------------------------------------------*/

Array<double> Reaction::GetMD()
{
    return MD;
}
